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CrossMark 


The combinatorial explosion of empirical parameters in tens of thousands presents a tremendous chal¬ 
lenge for extended Debye-Hiickel models to calculate activity coefficients of aqueous mixtures of the 
most important salts in chemistry. The explosion of parameters originates from the phenomenologi¬ 
cal extension of the Debye-Hiickel theory that does not take steric and correlation effects of ions and 
water into account. By contrast, the Poisson-Fermi theory developed in recent years treats ions and 
water molecules as nonuniform hard spheres of any size with interstitial voids and includes ion-water 
and ion-ion correlations. We present a Poisson-Fermi model and numerical methods for calculating 
the individual or mean activity coefficient of electrolyte solutions with any arbitrary number of ionic 
species in a large range of salt concentrations and temperatures. For each activity-concentration curve, 
we show that the Poisson-Fermi model requires only three unchanging parameters at most to well 
fit the corresponding experimental data. The three parameters are associated with the Born radius of 
the solvation energy of an ion in electrolyte solution that changes with salt concentrations in a highly 
nonlinear manner. Published by AIP Publishing. https://doi.Org/10.1063/l.5021508 


I. INTRODUCTION 

Thermodynamic modeling of aqueous electrolyte solu¬ 
tions plays an important role in chemical and biological sci¬ 
ences. 1_Ll Despite intense efforts in the past century, robust 
thermodynamic modeling of electrolyte solutions still presents 
a difficult challenge and remains a remote ambition in the 
extended Debye-Hiickel (DH) models due to the enormous 
number of parameters that need to be adjusted, carefully 
and often subjectively. 11 ' 1 ' For example, the Pitzer model 
requires 8 parameters for a ternary system and up to 8 tem¬ 
perature coefficients (parameters) for every Pitzer parameter 
in a temperature interval from 0 to about 200 °C. 11,13 It is 
indeed a frustrating despair (frustration on p. 11 in Ref. 9 
and despair on p. 301 in Ref. 1) that approximately 22 000 
parameters for combinatorial solutions of the most important 
28 cations and 16 anions in salt chemistry have to be extracted 
from the available experimental data for one temperature. 11 
The Pitzer model is still the most widely used DH model 
with unmatched precision for modeling aqueous electrolyte 
solutions over wide ranges of composition, temperature, and 
pressure. 13 

The Pitzer model and its variants 13 are all derived from 
the Debye-Hiickel theory 14 that in turn is based on a lin¬ 
ear Poisson-Boltzmann (PB) equation 11 although potentials 
calculated from PB near ions (for example) are often far 
beyond the linear range of the potential near ions or inter¬ 
faces. The PB equation treats ions as point charges with¬ 
out steric volumes and water molecules as a homogeneous 
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dielectric medium without steric volumes either and with a 
constant dielectric constant that neglects ion-water and ion-ion 
correlations. These simplifications give rise to the elegant, sim¬ 
ple, and useful DH theory. However, it is precisely because of 
the linearization and simplifications on steric and correlation 
effects that extended DH models have needed an explosion 
in the number of parameters in order to overcome the defi¬ 
ciencies (simplifications) of the classical Poisson-Boltzmann 
theory. The nonlinear PB equation was developed by Gouy and 
Chapman. 1516 

In the past few years, we have intensively investigated 
these two effects in a range of areas from electric dou¬ 
ble layers 1713 and ion activities 19 to biological ion chan¬ 
nels 18 - 20_24 and consequently developed an advanced theory— 
the Poisson-Fermi (PF) theory—that treats ions and water 
molecules as nonuniform hard spheres of any size with inter¬ 
stitial voids and includes many of the correlation effects of 
ions and water. We refer to our previous papers and refer¬ 
ences therein for a historical account of the literature of this 
theory. In Ref. 19, we proposed a PF model for calculating 
activity coefficients of individual ions in aqueous single NaCl 
and CaCL electrolyte solutions at the temperature 298.15 K. 
The model is further tested in this paper for eight 1:1 elec¬ 
trolytes (LiCl. LiBr, NaF, NaCl, NaBr, KF, KC1, and KBr), 
six 2:1 electrolytes (MgCL, Mg Bn, CaCL, CaBn, BaCL, 
and BaBro), one mixed electrolyte (NaCl + MgCL), one 1:1 
electrolyte (NaCl) at various temperatures from 298.15 to 
573.15 K, and one 2:1 electrolyte (MgCL) at various tempera¬ 
tures from 298.15 to 523.15 K, for which the experimental 
data were compiled by Valisko and Boda in Ref. 25 and 
Rowland et al. in Ref. 13 from various experimental sources in 
Refs. 26-35. 
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The PF model is developed to calculate individual ion 
activities for which experimental measurements and determi¬ 
nation, 10 " 16,37 interpretation of measurement data, 26,37-39 and 
comparison of different experimental methods 37,40 have been 
extensively investigated by Wilczek-Vera, Rodil, and Vera in 
the past two decades. PF results on mean activity coefficients 
can be compared with experimental measurements using the 
Debye-Hiickel equation of individual ion activities. 5 

In contrast to the Pitzer model, we show that all exper¬ 
imental data sets of individual or mean activity coefficients 
as a function of variable concentration in single electrolytes 
or mixtures at various temperatures can be well fitted by the 
PF model with only 3 parameters at most for each activity- 
concentration data curve. The model is characterized by three 
different domains, namely, the Bom ion, hydration shell, and 
remaining solvent domains in which the Bom ion domain is 
most crucial because all activities around an ion are mainly 
governed by the singular charge of the ion located at the cen¬ 
ter of the domain. The Born ion domain is defined by the Born 
radius of the solvated ion, which is unknown and changes with 
salt concentrations in a highly nonlinear manner. 

The three parameters characterize three orders of approx¬ 
imation of the Born radius in terms of ionic concentrations. 
Parameter 1 describes a correction of the experimental Born 
radius of a single ion in pure water without any other ions. 
Parameter 2 describes an adjustment of the unknown Born 
radius in electrolyte solution that accounts for the Debye 
screening effect, which is proportional to the square root of 
the ionic strength of the solution. Parameter 3 is an adjust¬ 
ment in the next order approximation beyond the DH treatment 
of ionic atmosphere. The physical origin of these param¬ 
eters is clear unlike that of most parameters in the Pitzer 
method. 11,41 It may even be possible in later work to calculate 
some of these parameters from more detailed versions of our 
model. 

Our approach to partition the free energy domain of a 
solvated ion into the above three sub-domains yields a better 
approximation to calculate the free energy since these sub- 
domains are determined by the experimental data of solvation 
and thus separate short- and long-range interactions of the ion 
in a more accurate way. This approach nevertheless incurs 
more complicated numerical methods for solving the nonlin¬ 
ear partial differential equations of the PF model in different 
domains with suitable interface conditions. 17 We therefore 
present numerical methods in detail for future verification and 
development of the present work. 

II. THEORY 

For an aqueous electrolyte solution with K species of ions, 
the Poisson-Fermi theory proposed in Refs. 18 and 21 treats 
all ions and water of any diameter as nonuniform hard spheres 
with interstitial voids between these spheres. The activity coef¬ 
ficient y,- of an ion of species i in the solution describes the 
deviation of the chemical potential of the ion from ideality (y, 
= 1). The excess chemical potential = k B T In y, can be 
calculated by 19,42 

pT = AG < - AG < = \<nm, AG - = l - qi 4>\ 0 ), ( 1 ) 



Solvent 

Domain 

£2 


FIG. 1. The model domain £2 is partitioned into the ion domain £2, (with radius 
R? on '), the hydration shell domain £2,/, (with radius RJ 1 '), and the remaining 
solvent domain Cl s . 


where kg is the Boltzmann constant, T is an absolute temper¬ 
ature, qt is the ionic charge of the hydrated ion (also denoted 
by ;), (Mr) is a potential function of spatial variable r in the 
domain Q = Q, U Q v /, U Q v shown in Fig. 1, Q, is the spherical 
domain occupied by the ion i, Q v /, is the hydration shell domain 
of the ion, Q s is the remaining solvent domain, 0 denotes the 
center (set to the origin) of the ion, <p( 0) is the value of (Mr) 
at r = 0, and (f>°( r) is a potential function when the solvent 
domain Q v does not contain any ion at all with pure water 
only. The potential function <p( r) can be found by solving the 
Poisson-Fermi equation 18 

(/2 V 2 - l) V • e(r)V0(r) - p(r), (2) 


I e s = e w eo in Qj, Ufi, [ 2 a, in £2 S /, U Q, 

e(r)= ' _ Jc = \ — 

[e,-= e ion e 0 in D,- (0 in Q, 


p(x) 


p s { r) = 2f =1 qkCk(r) in Q s 
0 in 


.Pi( r) = qMr - 0) in Q, 


Cjfc(r) = Cf exp \-/3 k <p(r) + ^S trc (r)j in O, 
5 ttc (r) = ln^J in D, 


(3) 

(4) 

(5) 

( 6 ) 


where ep is the vacuum permittivity, e w is the dielectric con¬ 
stant of bulk water, e lon is a dielectric constant in Q,, aj is the 
radius of a counterion of the ion ;, and <5(r - 0) is the delta 
function at the origin. 

The concentration function C, t(r) is described by a Fermi 
distribution (5), where C° is a constant bulk concentration for 
all k = 1, ..., K + 1, qx+i = 0, Pk = q k lk B T, v k = Ana^/3, 
vo = (2 k=i Vk ) /(^ + 1) an average volume of all kinds of 
hard spheres, .S’ uc (r) is called the steric potential, l B = I 
- Ef=/ v k Cf is a constant void fraction, T(r) = 1 - Tif=i 
v k C k (r) is a void fraction function, and K + 1 denotes water. 
The radii of Q, and the outer boundary of Q v /, are denoted by 
rBoi' 1 anc j /y/y respectively, whose values will be determined 
by experimental data. It is natural to choose the Born radius 
R /)V "" (not the ionic radius a,) as the radius of Q, 42 We consider 
both first and second shells of the ion. 41,44 
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The potential <p°( r) [in Eq. (1)] of the ideal system is 
obtained by setting p s (r) = 0 in (4), i.e., all particles in Q s do not 
electrostatically interact with each other since q k = 0 for all k. 
The domain Q. is chosen to be sufficiently large so that <p(r) = 0 
on the boundary of the domain dCl. The ideal potential </;°(r) 
is then a constant, i.e., AG° is a constant reference chemical 
potential independent of C®. 

The distribution (5) is of Fermi type since all concentration 
functions have an upper bound, i.e., C k (r) < l/vk for all particle 
species with any arbitrary (or even infinite) potential <f>( r) at any 
location r in the domain SI. 21 The Poisson-Fermi equation (2) 
and the Fermi distribution (5) reduce to the Poisson-Boltzmann 
equation and the Boltzmann distribution when l c = ,S' llc = 0, 
i.e., when the correlation and steric effects are not considered. 
The Boltzmann distribution C, t(r) = C® exp (-fii-(Mr)) would 
however diverge if <p(r ) tends to infinity. This is a major defi¬ 
ciency of the PB theory for modeling a system with strong local 
electric fields or interactions. 4 ^ If the correlation length l c + 0, 
the dielectric operator e = e v ( I - / r 2 V 2 J in Eq. (2) approximates 
the permittivity of the bulk solvent and the linear response of 
correlated ions 17 ' 20,46,47 and yields a dielectric function e(r) as 
an output of solving Eq. (2). 21 The exact value of e(r) at any 
r e O v /,UO v cannot he obtained from Eq. (2) hut can be approx¬ 
imated by the simple formula e(r) ~ e,-+C h, o(r)(e s -e,j/C ® 2Q 
since the water density function Ch^oM = CV+i (r) is an 
output of Eq. (5). This formula is only for visualizing (approxi¬ 
mately) the profile of e or e. It is not an input of calculation. The 
input is the correlation length l c in Eq. (3). 17 - 20 - 46 ' 47 The actual 
outputs are the numerical solutions of the partial differential 
equations and boundary conditions. 

The factor v k /vo multiplying the steric potential function 
A trc (r) in Eq. (5) is a modification of the unity used in our 
previous work. 19,21 The steric energy - —.S’ trc (r)k/j7’ 21 ’ 24 of a 
type k particle depends not only on the voidness (T(r)) (or 
equivalently crowding) at r but also on the volume Vk of the 
particle itself. If all Vk are equal (and thus v k = no), then all par¬ 
ticle species at any location r e fl s /, U Q s have the same steric 
energy, i.e., uniform particles are indistinguishable in steric 
energy. The steric potential is a mean-field approximation of 
Lennard-Jones (L-J) potentials that describe local variations 
of L-I distances (and thus empty voids) between any pair of 
particles. L-J potentials are highly oscillatory and extremely 
expensive and unstable to compute numerically. 21 Calcula¬ 
tions that involve L-J potentials or even truncated versions 
of L-J potentials must be extensively checked to be sure that 
results do not depend on irrelevant parameters. 

III. METHODS 

To avoid large errors in approximation caused by the 
delta function b(r - 0) in (4), the potential function can be 
decomposed as 17 - 48 - 49 

f0(r) + 0*(r) + 0 L (r)inO,- 

<P(r) = \ _ , (7) 

(0(r) in Oj/, U 

where <p*(r) = qtKAnei |r - 0 |) and </>(r) is found by solving 
(/ 2 V 2 - l) V • e s V0(r) = p(r) in U sh U Q s , (8) 

-V • e,V0(r) = 0 in O, (9) 


without the singular source term p,(r) = <r/,<5(r - 0 ) and with 
the interface conditions 

[0(r)l = 0 

_ for all r e dQ,-, 

[e(r)V0(r) • n] = e,V (0*(r) + 0 L (r))j • n 

(10) 

where n is an outward normal unit vector at r e dil, and the 
jump function [w(r)] = lim rj(i _, r u{r s i,) - lim ri _> r «(r,) with r s /, 
6 Q v /, and r, e Q,. 17 The potential function (f> L ( r) is the solution 
of the Laplace equation 

V 2 0 L (r) = 0 in Q, (11) 

with the boundary condition 

0 L (r) = <p*(r) on dQ,. (12) 

The evaluation of Green’s function (jr (r ) on dQ ., always yields 
finite numbers and thus avoids the singularity in the solution 
process. The desired solvation energy AG, in Eq. (1) (and thus 
the individual ionic activity coefficient y,) is then evaluated 

by 17 - 49 

AG, = k B T In y, = ±q { [^(0) + 0 L (O)] . (13) 

Since the interface <90, is a sphere centered at the origin, the 
Laplace potential <p L ( r) = qi/(4nejRf onl ) is a constant in O,, 
i.e., Eq. (11) has been exactly solved. 

The Poisson-Fermi equation (8) is a nonlinear fourth- 
order partial differential equation (PDE) in O s . Newton’s iter¬ 
ative method is usually used for solving nonlinear problems. 

We seek a sequence of approximate solutions {^ m (r)J by 
iteratively solving the linearized PF equation 

(/ 2 V 2 - l) V • eV^,„ - p'(0,„_i )(/>,„ 

= Ps(.<t>m- 1) - Ps(0m-l)0m-l in O s , (14) 

until a tolerable potential function </>m is reached, where 
0o(r) is a given initial guess potential function, p s (<p m -\) 
= 2f=t r), C?-'(r) = C® exp (- ^0 m -i(r) + |C-t 

(r^C-^r) = In(^), r m _!(r) = 1 - Zp} v k C™~Hr), 
p'Mm-i) = 2f=i i-Pkqk)C£- l (r), and p' s (</>) = p s (<p ). 

Note that the differentiation in p' s (<!>) is performed only with 
respect to (f >, whereas S tIC is treated as another indepen¬ 
dent variable although .S' trc depends on (f> as well. Therefore, 
p' s (<p) is not exact implying that this is an inexact Newton’s 
method. 50 

The fourth-order problem can be resolved by transforming 
Eq. (14) into two second-order PDEs 17 

e s (Z 2 V 2 - l) T'(r) = p(0,„_,) in a, h U Q„ (15) 

-e s V 2 (i>m(V) - p'((pm-\ )<pm{r) 

= -G'I'(r) - P\(p m - t)0m-t in kl sh U Q s (16) 

by introducing a density like variable T 1 = V 2 cf> for which the 
boundary condition is 17 


'P(r) = 0 on <9G V . 


(17) 
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Equations (9), (15), and (16) are coupled together in the entire 
domain Q with the jump conditions in (10). Note that linear 
PDEs (14)—(16) converge to the nonlinear PDE (8) if </>m con¬ 
verges to the exact solution / of Eq. (8) as M — » oo, i.e., the 
approximate potential 4>m( r) is sufficiently close to the exact 
potential Mr) for all r e O v ;, U Q s if the iteration number M is 
sufficiently large (M ~ 5-37 for this work with error tolerance 
10 " 3 ). 

The standard 7-point finite difference (FD) method is used 
to discretize all PDEs (9), (15), and (16), where the jump con¬ 
ditions in (10) are handled by the simplified matched interface 
and boundary (SMIB) method proposed in Ref. 17. For sim¬ 
plicity, the SMIB method is illustrated by the following ID 
linear Poisson equation (in x-axis): 


d 

dx 


d — 

e(x)—Mx) 

ax 


= fix) in Q 


( 18 ) 


with the jump condition 

\e(f>'] = -£j — /*(x) atx = £ = <90, n 0O s , (19) 

1 1 dx 

where Q = Q, U Q s , Q, = (0, /), Q s = (/, L),f(x) = 0 in O,, 
f(x) + 0 in and <p’ = fj^M x )- The corresponding cases to 
Eqs. (9), (15), and (16) in the y- and z-axis follow in a similar 
way. Let two FD grid points x/ and x/+i across the interface 
point / be such that x/ < / < x /+1 and / = (x/ + x/+i)/2 with 
Ax = x /+1 - xi = I A, a uniform mesh, for example, as used 
in this work. The FD equations of the SMIB method at x/ and 
x/+i are 


-0/-1 + (2 - C\)(pi - C20/+1 Co 

“ -S?- =/, + T?' <20) 

-di<pi + (2 - d2)(pi+\ - <f>i+2 , do 

'>- A? - =/M + A? (21) 


where 


2e, 


Cl 


e, + e. 


~,C2 


i 1 


-,c 0 


2e, 


-e,-Ax \cM\ 

£i + 

-e s Ax \eM\ 


d i = --—, <j?2 — ~---, d 0 = 

e, + e s ei + e s 6j + e. 


(jij is an approximation of <^(x/), and // =/(x/). Note that the 
jump value [ e( P'\ at £ is calculated exactly since the derivative 
of (p is given analytically. 

Since the steric potential takes particle volumes and voids 
into account, the shell volume V s h of the shell domain I//, can 
be determined by Eqs. (5) and (6) as 


^sh ~ — I * 1 

Vw 


D w \ 

—V- = ln 

VshC'J 


V sh -v w OJ\ 
VsX K J ’ 


( 22 ) 


where the occupancy (coordination) number OJ is given by 
experimental data. 43,44 The shell radius R S P' of Qis thus deter¬ 
mined. Note that the shell volume depends not only on OJ but 
also on the bulk void fraction T B , namely, on all salt and water 
concentrations (C 4! ). 

As discussed in Ref. 25, the solvation free energy of 
an ion i should vary with salt concentrations and can be 
expressed by a dielectric constant e(Cf) that depends on 


TABLE I. Values of model notations. 


Symbol 

Meaning 

Value 

Unit 

kjs 

Boltzmann constant 

1.38 x 10“ 23 

J/K 

T 

Temperature 

Table II 

K 

e 

Proton charge 

1.602 x 10“ 19 

C 

eo 

Permittivity of vacuum 

8.85 x 10“ 14 

F/cm 

^ion> 6 w 

Dielectric constants 

1, Table II 


lc — 2 <2/ 

Correlation length 

.. 

II 

Q 

CD 

O 

A 

Of 

In Eq. (22) 

jg43,44 


a Li + > a Na + > a K + 

Radii 

0.6, 0.95, 1.33 

A 

a Mg 2+ ’ fl Ca 2+ ’ fl Ba 2+ 

Radii 

0.65,0.99, 1.35 

A 

a F~ ’ a Cl~ ’ a Br ~» a H 2 0 

Radii 

1.36, 1.81, 1.95, 1.4 

A 

r>() p() p() 

Born radii in Eq. (24) 

1.3, 1.618, 1.95 

A 

pO pO pO 

Mg 2 +’ Ca 2 +’ Ba 2 + 

Born radii 

1.424, 1.708,2.03 

A 

p0 pO pO 
“f- ’ R cr' at- ■ 

Bom radii 

1.6,2.266,2.47 

A 


the bulk concentration C B of the ion. Therefore, the Born 

l 

energy 


AG Born 



8ne 0 R° 


(23) 


with the Born radius R^ in pure water should be modified 
with the concentration-dependent dielectric constant e(C B ). 
Equivalently, the Born radius in electrolyte solutions can be 
modified from by a simple formula 


R? orn {Cf) = 6(Cf)R°, 


0(Cf) = a[+a i 2 (cfj^ + a' (cf) 


/-R\3/2 


(24) 


where C, = C"/M is a dimensionless bulk concentration of 
type i ions, M is the molar concentration unit, and «', a' 2 , 
and crj are adjustable parameters for modifying the exper¬ 
imental Born radius to fit experimental activity coeffi¬ 
cients ji that change with the bulk concentration conditions 
C B of the ion. The Born radii R () in Table I are cited from 

l l 

Ref. 25, which are computed from the experimental hydration 
Helmholtz free energies of these ions given in Ref. 6. Numer¬ 
ical values in Tables I and II are all experimental data for 
which their values are kept fixed throughout calculations once 
chosen. 

The three parameters in Eq. (24) have physical or mathe¬ 
matical meaning unlike many parameters in the Pitzer model. 41 
Any model or numerical method incurs errors to approxi¬ 
mate a real system, i.e., it is impossible to obtain real Born 
radius R^ om (Cf) exactly. The first parameter a\ is an adjust¬ 
ment of the experimental Born radius R { . when C B = 0 
for all i. The second parameter a' 2 is an adjustment of 
Rf orn (Cf) that accounts for the real thickness of the ionic 
atmosphere (Debye length), which is proportional to the 
square root of the ionic strength V7 in the Debye-Hiickel the¬ 
ory. 3 The third parameter a/ is simply an adjustment in the 


TABLE II. Values of e w at various T. 51 


77K 

298.15 

373.15 

423.15 

473.15 

523.15 

573.15 


78.41 

55.51 

44.04 

38.23 

32.23 

25.07 
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next order approximation beyond the DH treatment of ionic 
atmosphere. 

We summarize the mathematical solution process for 
determining the activity of ionic solutions in the following 
algorithm. 

1. Solve Eqs. (9), (10), and (16) for <p with p' = T 7 = 0 (in 
pure water), R^ orn = R 9, and (f> L = ^ ( /(4^e,/?J ) ) to obtain 
AG° by Eq. (13) and then set (f >o = <p. 

2. Solve Eqs. (15) and (17) for ¥ with Rf orn in (24). 

3. Solve Eqs. (9), (10), and (16) for <f> m with (p L 
= qi/(4neiRf orn ) and then set (/>,„-\ = <p m . Go to 2 until 
convergence. 

4. Obtain the activity coefficient y, by Eq. (13). 

IV. RESULTS 

The PF results of ionic activity coefficients for eight 1:1 
electrolytes, six 2:1 electrolytes, one mixed electrolyte, one 1:1 
electrolyte at various temperatures, and one 2:1 electrolyte at 
various temperatures agree with the experimental data 2 ' 1-35 as 
shown in Figs. 2-6, respectively. The empirical parameters 
used to fit the experimental data are a j, a’ 2 , and a‘ 3 in Eq. (24), 
whose values are given in Table III from which we observe 
that the PF model requires only one to three parameters to fit 
those data. 

The mean activity coefficient yp OS Neg of a salt Pos^Neg^ 
is calculated via the formula In yp OS Neg = In yp os 
+ In y^eg , 5 where yp os and y^ eg are individual activity coef¬ 
ficients obtained by Eq. (13) for each i = Pos and Neg. For the 
mean activity coefficients of either ternary (Fig. 4) or binary 
(Figs. 5 and 6) systems, we only need to adjust 3 parameters 
of one cation (not all ions) as shown in Table III. 

The activity coefficients by the PF model are quite suc¬ 
cessful over a large range of temperatures and concentrations 
as shown in Figs. 4-6. We used the code of the density model 



FIG. 3. Individual activity coefficients of 2:1 electrolytes. Comparison of PF 
results with experimental data 26 on i = Pos 2+ (cation) and Neg - (anion) activity 
coefficients y,- in various [PosNeg 2 ] from 0 to 1.5M. 


developed by Mao and Duan 52 to convert the concentration 
unit from molality (mol Kg -1 ) to molarity (M = mol dm -3 ) 
by the standard formula as given in Ref. 52, where the den¬ 
sity model has been compared with thousands of measure¬ 
ments at high accuracy. The pressure values needed in the 
code at the corresponding temperatures were set to P = (a) 
1.01, (b) 1.01, (c) 15.48, (d) 39.59, and (e) 80.50 bars for 
Fig. 5 and (a) 1.01, (b) 1.01, (c) 4.73, and (d) 39.50 bars for 
Fig. 6. In Fig. 4, the ionic strength I = Cfzf and the ionic 
strength fraction y Mg ci 2 = 3m Mg ci 2 /(3m Mg ci 2 + «JNaCi) with 
/«M g c:i 2 and /M\ a ci being the molalities of MgCE and NaCl 



FIG. 2. Individual activity coefficients 
of 1:1 electrolytes. Comparison of PF 
results with experimental data 26 on i 
= Pos + (cation) and Neg - (anion) activ¬ 
ity coefficients y; in various [PosNeg] 
from 0 to 1.6M. 
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FIG. 4. Mean activity coefficients of mixed electrolytes. Comparison of PF 
results (curve) with experimental data (symbols) compiled in Ref. 13 (a) from 
Ref. 33 on mean activity coefficients y of NaCl as a function of the ionic 
strength (I) fraction yMgC^ °f MgCb in NaCl + MgCL mixtures at I = 6 mol 
Kg -1 and T = 298.15 K and (b) from Ref. 34 (circles) and Ref. 35 (squares) 
on y of NaCl as a function of the MgCl 2 molality in NaCl + MgCL> mixtures 
at [NaCl] = 6 mol Kg -1 and T = 298.15 K. 

in the mixture, respectively, where z, is the valence of type i 
ions. 

We observe from Table III that the approximate R? orn (Cf) 
(with salts) deviates from/?* 5 (without salts) only in the second 
to fourth decimal place, i.e., numerical values of y, are very 
sensitive to the decimal order of a' v a\, and a' 3 because the 
Born radius R?" ni (Cf) is very close to the origin 0 at which the 
singular charge in p,(r) = q,d(v - 0) is infinite. The approxi¬ 
mation of the shell radius R sh [or the coordination number ()' x 
in Eq. (22)], on the other hand, is much less significant than 
that of R? orn because the electric potential r/> PF (r) diminishes 
exponentially in the hydration shell region Q*/, as shown by 



FIG. 5. Mean activity coefficients of 1:1 electrolyte at various temperatures. 
Comparison of PF results (curves) with experimental data (symbols) compiled 
in Ref. 13 from Refs. 27-29 on mean activity coefficients y of NaCl in [NaCl] 
from 0 to 6 mol Kg -1 at T = (a) 298.15, (b) 373.15, (c) 473.15, (d) 523.15, 
(e) 573.15 K. 



FIG. 6. Mean activity coefficients of 2:1 electrolyte at various temperatures. 
Comparison of PF results (curves) with experimental data (symbols) compiled 
in Ref. 13 from Refs. 30-32 on mean activity coefficients y of MgCl 2 in 
[MgCl 2 ] from 0 to 6 mol Kg -1 at T = (a) 298.15, (b) 373.15, (c) 423.15, (d) 
523.15 K. 

the profile of <^ PF (r) in Fig. 7. The values of a‘ v a' 2 , and a‘ 3 
for each activity-concentration curve were obtained by first 
tuning three values of G(Cf) in Eq. (24) to match three data 

points (yC?, lnyy) with three different concentrations C?,y 
= 1, 2, 3, and then solving the three unknowns a' { , a' 2 , and a' 3 
using three known 0(C?) values. For example, for the / = Li + 
curve in Fig. 2(a), the selected experimental data points are 

( A /C?,ln yu) = (0.315, -0.192), (1, -0.007), and (1.577, 0.57) 

y l J J 

and the corresponding tuned 6{Cf-) are 0.9996, 1.0013, and 
1.0043. 

The PF model can provide more physical details near 
the solvated ion (Ca 2+ , for example) in a strong electrolyte 
([CaCy = 2M) such as (1) the dielectric function e(r) with 
its varying permittivity, (2) variable water density CH 2 o( r ), 
(3) concentration of counterion C cr (r), (4) electric potential 
0 PF (r), and (5) the steric potential ,S' trc (r) all shown in Fig. 7. 
The steric potential is small because the configuration of par¬ 
ticles (voids between particles) does not vary too much from 
the solvated region to the bulk region. Nevertheless, it has sig¬ 
nificant effect on the variation of mean-field water densities 
CH 2 o( r ) and hence on the dielectric function F(r) in the hydra¬ 
tion region. Note that e(r) is an output, not an input of the 
model. 

The strong electric potential 0 PF (r) in the Born cav¬ 
ity O i (with R? orn (Cf) = 1.7130 A) and the water density 
CH 2 o( r ) in the hydration shell Q v /, (with R' h 2+ = 5.0769 A) 
are the most important factors allowing the PF results to 
match the experimental data. The ion and shell domains 
are the crucial region to study ion activities. For exam¬ 
ple, Fraenkel’s theory is entirely based on this region—the 
so-called smaller-ion shell region. 41 The steric energy of 
water molecules modified by the factor Vk+i/vo in Eq. (5) 
leads to significant changes of ChicM and <T(r) profiles in 
Fig. 7 as compared with those in Fig. 5 in our previous 
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paper. 
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TABLE III. Values of a' { , a' 2 , and a' 2 in Eq. (24). 


Figures 

i 

a\ 

ry 1 

2 

ry 1 

J 

Figures 

i 

a\ 

a 2 

ry 1 

J 

2(a) 

Li+ 

0.99913 

0.00069 

0.00009 

3(c) 

Ca 2+ 

0.998 86 

0.00046 

0.00011 

2(a) 

Cl - 

0.998 93 

-0.00008 


3(c) 

cr 

0.998 77 

-0.00060 

0.000 12 

2(b) 

Li+ 

0.999 58 

-0.000 19 

0.000 15 

3(d) 

Ca 2+ 

0.998 86 

0.00099 

0.000 17 

2(b) 

Br - 

0.998 22 

0.001 07 


3(d) 

Br - 

0.999 20 

-0.001 98 

0.000 16 

2(c) 

Na+ 

0.999 10 



3(e) 

Ba 2+ 

0.998 44 

0.00011 

0.000 10 

2(c) 

F - 

0.99933 

-0.00029 


3(e) 

cr 

0.998 87 

-0.00058 

0.00001 

2(d) 

Na + 

0.999 27 

0.00026 

0.00004 

3(f) 

Ba 2+ 

0.998 51 

0.00054 

0.00008 

2(d) 

Cl - 

0.998 40 



3(f) 

Br - 

0.999 26 

-0.001 45 

0.00018 

2(e) 

Na + 

0.999 62 

-0.000 38 

0.000 10 

4(a) 

Na + 

1.005 81 

-0.000 13 


2(e) 

Br - 

0.998 70 

-0.000 17 

0.00004 

4(b) 

Na + 

1.005 27 

0.00042 

0.000 19 

2(f) 

K + 

0.999 34 

-0.001 20 

0.00007 

5(a) 

Na + 

0.998 1 


0.000 1 

2(f) 

F - 

0.99904 

0.00013 

0.00004 

5(b) 

Na + 

0.997 1 

0.000 3 

0.000 1 

2(g) 

K + 

0.999 29 

-0.001 22 

0.00004 

5(c) 

Na + 

0.9945 

-0.0007 

0.000 1 

2(g) 

Cl - 

0.998 97 

-0.000 12 

0.00003 

5(d) 

Na + 

0.992 5 

-0.002 8 

0.000 1 

2(h) 

K + 

0.999 31 

0.000 13 


5(e) 

Na + 

0.987 0 

-0.0042 

0.001 0 

2(h) 

Br - 

0.99945 

-0.001 75 

-0.00006 

6(a) 

Mg 2+ 

0.998 8 

0.0002 

0.0002 

3(a) 

Mg 2+ 

0.99918 

0.00044 

0.00011 

6(b) 

Mg 2+ 

0.998 9 

-0.0004 

0.000 3 

3(a) 

Cl - 

0.998 93 

-0.00051 

0.000 10 

6(c) 

Mg 2+ 

0.998 3 

-0.001 4 

0.0005 

3(b) 

Mg 2+ 

0.999 10 

0.00063 

0.000 15 

6(d) 

Mg 2+ 

0.996 1 

-0.002 0 

0.000 3 

3(b) 

Br - 

0.998 88 

-0.00065 

0.000 18 






Default values: orj 

= l,«j = 

0, and CTj = 0. 









FIG. 7. Dielectric function e(r) (denoted by s in the fig¬ 
ure), water density Cn ,o(r) (Ch,o), Cl - concentration 
Cq-( r) ([Cl - ]), electric potential <b n (r) ( tp ), and steric 
potential .V trc (r) (.S’ Irc ) profiles near the solvated ion Ca 2+ 
at [CaCy = 2M, where r is the distance from the center 
of Ca 2+ in angstrom. 


V. CONCLUSION 

A Poisson-Fermi model for calculating activity coeffi¬ 
cients of aqueous single or mixed electrolyte solutions in a 
large range of concentrations and temperatures has been pre¬ 
sented and tested by a set of experimental data. The model was 
shown to well fit experimental data with only three adjustable 
parameters at most for each activity-concentration curve. 
The adjustable parameters correspond to different orders of 
approximation of the unknown Born radius of solvation energy 
that depends on salt concentrations in a highly complex and 
nonlinear way. Nevertheless, the values of these parameters 
have been shown to deviate slightly in decimal digits from 
that of the experimental Born radius in pure water. These 
parameters are physically explained and can be easily veri¬ 
fied in future studies for the same or different solutions of the 
present work. The model requires very few parameters because 
it is based on an advanced continuum theory that accounts 
for steric and correlation effects of ions and water with 


interstitial voids between nonuniform hard spheres. It also 
deals with short- and long-range interactions by partitioning 
the model domain into the ion, hydration shell, and the remain¬ 
ing solvent sub-domains. Numerical methods were also given 
to show how to solve different equations on different sub- 
domains that describe different physical properties of an ion 
in electrolyte solutions. 
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